n=2.05;

eqX=m x''[t]==-G M m /Sqrt[x[t]^2+y[t]^2+z[t]^2]^(n+1) x[t]

eqY=m y''[t]==-G M m /Sqrt[x[t]^2+y[t]^2+z[t]^2]^(n+1) y[t]

eqZ=m z''[t]==-G M m /Sqrt[x[t]^2+y[t]^2+z[t]^2]^(n+1) z[t]

G=1

m=1 (* Masa próbna, skraca się ! *)

M=1

d=1

ω=Sqrt[G M/d^3]

v0=ω d/3

E0=-G M m/d+1/2 m v0^2 (* Energia z warunków początkowych, zachowana podczas ruchu *)

a = -G m M/2/E0 (* Wzór na wielką półoś elipsy, słuszny dla potencjału 1/r *)

T=2Pi/Sqrt[G M/a^3]

j=m v0 d

e=Sqrt[1+2 E0 j^2/G^2/m^3/M^2] (* Mimośród elipsy, słuszny w potencjale 1/r *)

rozw=NDSolve[
{eqX,eqY,eqZ,
x[0]==d,y[0]==0,
x'[0]==0,y'[0]==-v0,
z[0]==0,z'[0]==0},
{x,y,z},{t,0,8T},
Method->"BDF"
][[1]]

X=x/.rozw;
Y=y/.rozw;
Z = z/.rozw;

X[1]

8.0T

Plot[Z[t],{t,0,8T}]

ParametricPlot[{X[t],Y[t]},{t,0,8T},
PlotRange->1.2d{{-1,1},{-1,1}},AxesOrigin->{0,0,0},PlotStyle->{Black}(*,PlotLabel->"1/"r^n*),Axes->False,BaseStyle->{FontSize->24},PlotPoints->128,
Prolog->Point[{0,0}],
Epilog->{Arrow[BSplineCurve[Table[1.2 d{X[τ],Y[τ]},{τ,0,8T,T}]]]}(*,
Epilog->{Arrow[Table[1.2 d{X[τ],Y[τ]},{τ,0,8T,T}]]}*)
]

Na ,,oko’’  widać, że elipsa zachowuje swój rozmiar (wielką oś 2a) oraz mimośród. Natomiast obraca się. Zjawisko to określamy precesją peryhelium i jest to typowy przypadek zaburzenia, które zachowuje eliptyczność orbity. Natomiast jej parametry zmieniają się powoli. Rysunek pokazuje zjawisko z ogromną przesadą, np. słynna precesja peryhelium Merkurego to 5600’’ (sekund kątowych) na 100 lat, w tym zaledwie 43’’ efektów OTW. Wizualnie oczywistej precesji można spodziewać się jedynie w rejonach czarnych dziur i gwiazd neutronowych. 

Pomysł na projekt: czy powyższe stwierdzenie jest prawdziwe, czy wynika z GR-centrycznych poglądów? Co np. z planetą orbitującą konaktowy układ podwójny gwiazd, albo satelitą na orbicie planetoidy o dziwacznym kształcie. S. Lem w Dziennikach gwiazdowych opisywał  sytuację, w której bohater zgubił wołowinę i klucz francuski. Następnie z nudów obliczył ich elementy orbitalne i wyszło mu że, przez 6 milionów lat wołowina będzie wyprzedzała klucz...

Quantity[5600.0 ,"arcseconds"]

UnitConvert[%,"Degrees"]

Bardziej naukowe podejście niż ocena precesji na oko, polega na użyciu zachowanych w potencjale 1/r wektorów momentu pędu J oraz Rungego-Lenza A. W ogólności nie są one zachowane, ale na skali jednej orbity moża je uśrednić i potraktować jako wolnozmienne funkcje.

RungeLenz=Cross[ m {X'[t],Y'[t],Z'[t]},
Cross[{X[t],Y[t],Z[t]},m{X'[t],Y'[t],Z'[t]}]]-G m^2 M{X[t],Y[t],Z[t]}/Norm[{X[t],Y[t],Z[t]}] ;

eRL=RungeLenz/G m^2 M;

e/Norm[eRL]-1/.t->T (* Sprawdzenie zgodności mimośrodu e z war. początkowych i z wektora Rungego-Lenza *)

Plot[Norm[eRL],{t,0,8T},PlotRange->{0,1}]

Plot[VectorAngle[eRL,{1,0,0}],{t,0,8T},PlotRange->{-Pi,Pi}]

J=m Cross[{X[t],Y[t],Z[t]},{X'[t],Y'[t],Z'[t]}];

Plot[J[[3]],{t,0,8T},PlotRange->All,PlotStyle->{Red,Green,Blue}]

Plot[RungeLenz,{t,0,8T},PlotRange->All,PlotStyle->{Red,Green,Blue}]

Plot[RungeLenz.J,{t,0,8T}]